/**
 * @file   Element.h
 * @author ubuntu <dadi@ubuntu>
 * @date   Mon Jun 14 21:25:14 2021
 * 
 * @brief  
 * 
 * 
 */

#ifndef __DADI_ELEMENT__
#define __DADI_ELEMENT__

#include "Point.h"
#include "Matrix.h"
#include <valarray>
#include <iostream>
#include <cstdlib>
#include <fstream>

#define TEMPLATE template <int DIM, typename T>
#define TEMPLATE2D template <typename T>


TEMPLATE class Element;

TEMPLATE2D
class TriElement;

TEMPLATE class Element
{
protected:
    std::valarray<Point<DIM, T> > local_dofs;
    std::valarray<Point<DIM, T> > global_dofs;
    std::valarray<Point<DIM, T> > quad_points;
    std::valarray<T> quad_weights;
    int algebraic_accuracy;
public:
    Element() {algebraic_accuracy=0;};
    void set_n_dofs(size_t _n);
    size_t get_n_dofs() const;
    void set_local_dof(size_t _i, const Point<DIM, T> &_pnt);
    void set_global_dof(size_t _l, const Point<DIM, T> &_pnt);
    const Point<DIM, T>& get_local_dof(size_t _i) const;
    const Point<DIM, T>& get_global_dof(size_t _i) const;
    const std::valarray<Point<DIM, T> >& get_local_dofs() const;
    const std::valarray<Point<DIM, T> >& get_global_dofs() const;
    size_t get_n_quad_points();
    const std::valarray<Point<DIM, T> >& get_quad_points() const;
    const std::valarray<T>& get_quad_weights() const;
    size_t get_algebraic_accuracy() const;
    virtual T lambda_function(size_t _i, const Point<DIM, T> &_pnt) const = 0;
    virtual std::valarray<T> lambda_gradient(size_t _i, const Point<DIM, T> &_pnt) const = 0;
    virtual Point<DIM, T> local2global(const Point<DIM, T> &_pnt) const = 0;
    virtual T basis_function(size_t _i, const Point<DIM, T> &_pnt) const = 0;
    virtual std::valarray<T> basis_gradient(size_t _i, const Point<DIM, T> &_pnt) const = 0;
};

TEMPLATE
void Element<DIM, T>::set_n_dofs(size_t _n)
{
    local_dofs.resize(_n);
    global_dofs.resize(_n);
};

TEMPLATE
size_t Element<DIM, T>::get_n_dofs() const
{
    size_t n_local_dofs = local_dofs.size();
    size_t n_global_dofs = global_dofs.size();
    if (n_local_dofs != n_global_dofs)
    {
	std::cerr << "No. of dofs error." << std::endl;
	std::exit(-1);
    }
    return n_local_dofs;
};

TEMPLATE
size_t Element<DIM, T>::get_n_quad_points()
{
    size_t n_q_pnts = quad_points.size();
    size_t n_q_wei = quad_weights.size();
    if (n_q_wei != n_q_pnts)
    {
	std::cerr << "No. of quad points error." << std::endl;
	std::exit(-1);
    }
    return n_q_pnts;
};

TEMPLATE
size_t Element<DIM, T>::get_algebraic_accuracy() const
{
    return algebraic_accuracy;
};


TEMPLATE
void Element<DIM, T>::set_global_dof(size_t _i, const Point<DIM, T> &_pnt)
{
    global_dofs[_i] = _pnt;
    global_dofs[_i].set_index(_pnt.get_index());
};

TEMPLATE
const Point<DIM, T>& Element<DIM, T>::get_local_dof(size_t _i) const
{
    return local_dofs[_i];
};

TEMPLATE
const Point<DIM, T>& Element<DIM, T>::get_global_dof(size_t _i) const
{
    return global_dofs[_i];
};


TEMPLATE
const std::valarray<Point<DIM, T> >& Element<DIM, T>::get_local_dofs() const
{
    return local_dofs;
};

TEMPLATE
const std::valarray<Point<DIM, T> >& Element<DIM, T>::get_global_dofs() const
{
    return global_dofs;
};

TEMPLATE
const std::valarray<Point<DIM, T> >& Element<DIM, T>::get_quad_points() const
{
    return quad_points;
};

TEMPLATE
const std::valarray<T>& Element<DIM, T>::get_quad_weights() const
{
    return quad_weights;
};

TEMPLATE
void Element<DIM, T>::set_local_dof(size_t _i, const Point<DIM, T> &_pnt)
{
    local_dofs[_i] = _pnt;
    local_dofs[_i].set_index(_i);
};

TEMPLATE2D
class TriElement : public Element<2, T>   
{
private:
    size_t n_dofs = 3;
public:
    TriElement();
    TriElement(int _n_dofs);
    void read_quad_info(const char* filename, size_t _acc);
    T lambda_function(size_t _i, const Point<2, T> &_pnt) const;
    std::valarray<T>  lambda_gradient(size_t _i, const Point<2, T> &_pnt) const;
    Point<2, T> local2global(const Point<2, T> &_pnt) const;
    Point<2, T> global2local(const Point<2, T> &_pnt) const;
    T global2local_Jacobi_det(const Point<2, T> &_pnt) const;
    Matrix<T> global2local_Jacobi(const Point<2, T> &_pnt) const;
    Matrix<T> local2global_Jacobi(const Point<2, T> &_pnt) const;
    T get_volume() const;
};

TEMPLATE2D
TriElement<T>::TriElement() : Element<2, T>()
{
    Element<2, T>::set_n_dofs(n_dofs);
    Point<2, T> pnt({0.0, 0.0});
    Element<2, T>::set_local_dof(0, pnt);
    pnt = {1.0, 0.0};
    Element<2, T>::set_local_dof(1, pnt);
    pnt = {0.0, 1.0};
    Element<2, T>::set_local_dof(2, pnt);
};

TEMPLATE2D
TriElement<T>::TriElement(int _n_dofs) : Element<2, T>()
{
    if (_n_dofs < 3)
	n_dofs = 3;
    else
	n_dofs = _n_dofs;
    Element<2, T>::set_n_dofs(n_dofs);
    Point<2, T> pnt({0.0, 0.0});
    Element<2, T>::set_local_dof(0, pnt);
    pnt = {1.0, 0.0};
    Element<2, T>::set_local_dof(1, pnt);
    pnt = {0.0, 1.0};
    Element<2, T>::set_local_dof(2, pnt);
};

TEMPLATE2D
T TriElement<T>::lambda_function(size_t _i, const Point<2, T> &_pnt) const
{
    double xi = _pnt[0];
    double eta = _pnt[1];
    switch(_i)
    {
    case 0:
	return (1 - xi - eta);
    case 1:
	return xi;
    case 2:
	return eta;
    }
};

TEMPLATE2D
std::valarray<T> TriElement<T>::lambda_gradient(size_t _i, const Point<2, T> &_pnt) const
{
    switch(_i)
    {
    case 0:
	return std::valarray<T>({-1.0, -1.0});
    case 1:
	return std::valarray<T>({1.0, 0.0});
    case 2:
	return std::valarray<T>({0.0, 1.0});;
    }
};

TEMPLATE2D
Point<2, T> TriElement<T>::local2global(const Point<2, T> &_pnt) const
{
    Point<2, T> global_pnt({0.0, 0.0});
    const std::valarray<Point<2, T> > &gdofs = Element<2, T>::get_global_dofs();
    for (size_t i = 0; i < 3; i++)
	global_pnt += gdofs[i] * lambda_function(i, _pnt); 
    return global_pnt;
};    

TEMPLATE2D
Point<2, T> TriElement<T>::global2local(const Point<2, T> &_pnt) const
{
    std::valarray<T> global({_pnt[0]-(this->global_dofs)[0][0],_pnt[1]-(this->global_dofs)[0][1]});
    std::valarray<T> local=local2global_Jacobi(_pnt).transpose()*(global);
    Point<2,T> local_pnt({local[0],local[1]});
    return local_pnt;
};

TEMPLATE2D
Matrix<T> TriElement<T>::global2local_Jacobi(const Point<2, T> &_pnt) const //d(x,y)/d(u,v)
{
//    const Point<2, T> &pnt0 = Element<2, T>::degree_of_freedoms[0][1];
//    const Point<2, T> &pnt1 = Element<2, T>::degree_of_freedoms[1][1];
//    const Point<2, T> &pnt2 = Element<2, T>::degree_of_freedoms[2][1];
    const std::valarray<Point<2, T> > &gdofs = Element<2, T>::get_global_dofs();
    Matrix<T> result(2, 2);
    result(0, 0) = gdofs[1][0] - gdofs[0][0];
    result(0, 1) = gdofs[1][1] - gdofs[0][1];
    result(1, 0) = gdofs[2][0] - gdofs[0][0];
    result(1, 1) = gdofs[2][1] - gdofs[0][1];
    return result;
};

TEMPLATE2D
T TriElement<T>::global2local_Jacobi_det(const Point<2, T> &_pnt) const
{
    // const Point<2, T> &pnt0 = Element<2, T>::degree_of_freedoms[0][1];
    // const Point<2, T> &pnt1 = Element<2, T>::degree_of_freedoms[1][1];
    // const Point<2, T> &pnt2 = Element<2, T>::degree_of_freedoms[2][1];
    // return ((pnt1[0] - pnt0[0]) * (pnt2[1] - pnt0[1])
    // 	    - (pnt1[1] - pnt0[1]) * (pnt2[0] - pnt0[0]));
    Matrix<T> J = global2local_Jacobi(_pnt);
    return J(0, 0) * J (1, 1) - J(0, 1) * J(1, 0);
};

TEMPLATE2D
T TriElement<T>::get_volume() const
{
    const Point<2, T> &pnt0 = Element<2, T>::local_dofs[0];
    const Point<2, T> &pnt1 = Element<2, T>::local_dofs[1];
    const Point<2, T> &pnt2 = Element<2, T>::local_dofs[2];
    return 0.5 * ((pnt1[0] - pnt0[0]) * (pnt2[1] - pnt0[1])
		  - (pnt1[1] - pnt0[1]) * (pnt2[0] - pnt0[0]));
    
};


TEMPLATE2D
Matrix<T> TriElement<T>::local2global_Jacobi(const Point<2, T> &_pnt) const //d(u,v)/d(x,y)
{
    Matrix<T> result(2, 2);
    Matrix<T> J = global2local_Jacobi(_pnt);
    T det = J(0, 0) * J (1, 1) - J(0, 1) * J(1, 0); 
    result(0, 0) = J(1, 1);
    result(0, 1) = -J(0, 1);
    result(1, 0) = -J(1, 0);
    result(1, 1) = J(0, 0);
    result /= det;
    return result;
};    

TEMPLATE2D
void TriElement<T>::read_quad_info(const char* filename, size_t _acc)
{
    std::cout << "Reading quadrature data ..." << std::endl;	
    std::ifstream is(filename);  
    char text[64];
    size_t n_tbl, n_qpoints, acc, next_acc;
    is >> n_tbl;
    if (n_tbl < _acc)
    {
	std::cerr << "The required algebraic accuracy " << _acc
		  << " is too high,"
		  << " and the actual algebraic accuracy is "
		  << n_tbl << std::endl;;
	acc = n_tbl;
    }
    else
	acc = _acc;
    do
    {
	is >> next_acc;
	if (next_acc >= acc)
	{
	    is >> n_qpoints;
	    Element<2, T>::quad_points.resize(n_qpoints);
	    Element<2, T>::quad_weights.resize(n_qpoints);
	    for (size_t i = 0; i < n_qpoints; i++)
	    {
		T x, y, w;
		is >> x >> y >> w;
		Element<2, T>::quad_points[i][0] = x;
		Element<2, T>::quad_points[i][1] = y;
		Element<2, T>::quad_weights[i] = w;
	    }
	    Element<2, T>::algebraic_accuracy = acc;
	    break;
	}
	else
	{
	    is >> n_qpoints;
	    for (size_t i = 0; i < n_qpoints; i++)
		is.getline(text, 64);
	    is.getline(text, 64);
	}
    } while (1);
    std::cout << "OK!" << std::endl;	
   
};

TEMPLATE2D
class P1Element : public TriElement<T>   
{
public:
    P1Element(){};
    T basis_function(size_t _i, const Point<2, T> &_pnt) const;
    std::valarray<T> basis_gradient(size_t _i, const Point<2, T> &_pnt) const;
};

TEMPLATE2D
T P1Element<T>::basis_function(size_t _i, const Point<2, T> &_pnt) const
{
    double xi = _pnt[0];
    double eta = _pnt[1];
    switch(_i)
    {
    case 0:
	return (1 - xi - eta);
    case 1:
	return xi;
    case 2:
	return eta;
    }
};

TEMPLATE2D
std::valarray<T> P1Element<T>::basis_gradient(size_t _i, const Point<2, T> &_pnt) const
{
    switch(_i)
    {
    case 0:
	return std::valarray<T>({-1.0, -1.0});
    case 1:
	return std::valarray<T>({1.0, 0.0});
    case 2:
	return std::valarray<T>({0.0, 1.0});;
    }
};

TEMPLATE2D
class P2Element : public TriElement<T>   
{
public:
    P2Element();
    T basis_function(size_t _i, const Point<2, T> &_pnt) const;
    std::valarray<T> basis_gradient(size_t _i, const Point<2, T> &_pnt) const;
};

TEMPLATE2D
P2Element<T>::P2Element() : TriElement<T> (6)
{
    Point<2, T> pnt({0.5, 0.5});
    Element<2, T>::set_local_dof(3, pnt);
    pnt = {0.0, 0.5};
    Element<2, T>::set_local_dof(4, pnt);
    pnt = {0.5, 0.0};
    Element<2, T>::set_local_dof(5, pnt);    
};

TEMPLATE2D
T P2Element<T>::basis_function(size_t _i, const Point<2, T> &_pnt) const
{
    double xi = _pnt[0];
    double eta = _pnt[1];
    switch(_i)
    {
    case 0:
	return (1 - xi - eta) * (1.0 - 2.0 * xi - 2.0 * eta);
    case 1:
	return xi * (2.0 * xi - 1.0);
    case 2:
	return eta * (2.0 * eta - 1.0);
    case 3:
	return 4.0 * xi * eta;
    case 4:
	return 4.0 * (1.0 - xi - eta) * eta;
    case 5:
	return 4.0 * (1.0 - xi - eta) * xi;
    }
};

TEMPLATE2D
std::valarray<T> P2Element<T>::basis_gradient(size_t _i, const Point<2, T> &_pnt) const
{
    double xi = _pnt[0];
    double eta = _pnt[1];
    double phi_xi, phi_eta;

    switch(_i)
    {
    case 0:
	phi_xi = phi_eta = -3.0 + 4.0 * xi + 4.0 * eta;
	return std::valarray<T>({phi_xi, phi_eta});
    case 1:
	phi_xi = 4.0 * xi - 1.0;
	phi_eta = 0.0;
	return std::valarray<T>({phi_xi, phi_eta});
    case 2:
	phi_xi = 0.0;
	phi_eta = 4.0 * eta - 1.0;
	return std::valarray<T>({phi_xi, phi_eta});;
    case 3:
	phi_xi = 4.0 * eta;
	phi_eta = 4.0 * xi;
	return std::valarray<T>({phi_xi, phi_eta});;
    case 4:
	phi_xi = -4.0 * eta;
	phi_eta = 4.0 - 4.0 * xi - 8.0 * eta;
	return std::valarray<T>({phi_xi, phi_eta});;
    case 5:
	phi_xi = 4.0 - 4.0 * eta - 8.0 * xi;
	phi_eta = -4.0 * xi;
	return std::valarray<T>({phi_xi, phi_eta});;
    }
};

#undef TEMPLATE
#undef TEMPLATE2D

#else
//DO NOTHING.
#endif
